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■ In this paper it has been described how to use the unitary dynamics 
O ! of quantum mechanics to solve the prime factorization problem on a spin 
^ I ensemble without any quantum entanglement. The ensemble quantum com- 
putation for the prime factorization is based on the basic principle that both 
a closed quantum system and its ensemble obey the same unitary dynamics 

■ of quantum mechanics if there is not any decoherence effect in both the quan- 
^ , tum system and its ensemble. It uses the NMR multiple-quantum measure- 
Q [ ment techniques to output the quantum computational results that are the 
2 I inphase multiple-quantum spectra of the spin ensemble. It has been shown 

■ that the inphase NMR multiple-quantum spectral intensities used to search 
p ■ for the period of the modular exponential function may reduce merely in a 

. polynomial form as the qubit number of the spin ensemble. The time evolu- 

^ I tion process of the modular exponential operation on the quantum computer 

^ I obeys the unitary dynamics of quantum mechanics and hence the computa- 

CT' tional output is governed by the quantum dynamics. This essential difference 

^ ■ between the quantum computer and the classical one could be the key point 

^ . for the quantum computation outperforming the classical one in the prime 

factorization on a spin ensemble without any quantum entanglement. It has 
been shown that the prime factorization based on the quantum dynamics on 
a spin ensemble is locally efficient at least. This supports the conjecture that 
the quantum dynamics could play an important role for the origin of power 
of quantum computation and quantum entanglement could not be a unique 
resource to achieve power of quantum computation in the prime factoriza- 
tion. 
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1. Introduction 

The prime factorization is an important problem that has made it a rapid 
development for the quantum computation and quantum information science. 
The Shor's quantum algorithm [1, 2] proposed first in 1994 to factorize effi- 
ciently a large composite integer can provide a possibility to break down the 
current public key cryptography such as the RSA cryptosystem. This fact 
has stimulated a great interesting in the quantum computation and informa- 
tion science and has promoted a large advance in the quantum computation 
[3, 4, 5, 6, 7]. This factorization algorithm could not come along for a short 
time most because the present-day quantum systems have not an enough 
long decoherence time to run the algorithm, although a preliminary exper- 
iment verification for the Shor's algorithm on an NMR quantum computer 
was reported [8]. The factoring algorithm is based on a pure-state quantum 
system. It has been suggested due to this powerful quantum algorithm that 
the exponential speedup power of the quantum computation over the classi- 
cal counterpart could be attributed to quantum entanglement of a quantum 
system [9]. One reason for it is that quantum entanglement is a uniquely 
feature differing the quantum effect from the classical effect and most power- 
ful quantum algorithms nowadays involve in the quantum entanglement [9]. 
Another is the well-known fact that quantum entanglement plays a key im- 
portant role in quantum communication [10, 11]. However, it has never been 
proved rigorously that quantum entanglement is the sole origin of power of 
quantum computation, and a number of recent works [12, 13] have showed 
that the power of quantum computation may not originate from quantum 
entanglement, although no work shows so far that the exponential speedup 
in the factoring algorithm may be independent of quantum entanglement. 
Very recently, an improved factoring algorithm has been proposed [14]. It 
has been shown that a quantum system consisting of an auxiliary pure-state 
qubit and log2 N mixed qubits is still sufficient to implement efficiently the 
prime factorization [14, 15]. But it also has been argued [14, 16] that quantum 
entanglement could play an important role in achievement of the exponen- 
tial speedup in the algorithm since there still exists quantum entanglement 
in such a system. The exponential speedup achieved on such a system really 
does not provide any certain answer whether or not quantum entanglement 
plays an important role in the exponential speedup of quantum computation 
over the classical computation. 

Fortunately, there are a lot of quantum ensembles in nature in which there 



2 



is not any quantum entanglement. These quantum ensembles include the 
conventional NMR nuclear spin ensembles at room temperature [17], which 
are also macroscopic quantum ensembles [18]. Quantum entanglement in a 
spin ensemble may be controlled by temperature of the spin ensembles. It is 
easy to keep any mixed state of a spin ensemble even with a larger number of 
qubits in a nonentanglement state by setting the spin ensemble at a higher 
temperature [17], but temperature of a spin ensemble should be as low as 
possible in order to make the NMR signal-to-noise ratio high enough for any 
NMR experiments. Such an ensemble without any quantum entanglement 
could be a typical system to judge whether or not quantum entanglement 
is the origin of power of quantum computation. Actually, if any quantum 
algorithm such as the factoring algorithm could be implemented efficiently 
in such a quantum ensemble without any quantum entanglement one could 
conclude certainly that quantum entanglement is not the unique origin of 
power of quantum computation. This is one of the reasons why the NMR 
spin ensembles are chosen as the typical systems to study the origin of power 
of quantum computation in the paper. Another reason is that the NMR spin 
ensembles usually have a long relaxation time and are simple and easy to 
be controlled and manipulated at will in experiments, and there are a large 
number of well-developed experimental techniques in the NMR spectroscopy 
[19, 20] which all can be adopted in NMR quantum computation. 

The unitary dynamic method of quantum mechanics has been proposed 
to solve efficiently the quantum search problem and the hard NP-problems 
[21, 22, 23, 24]. It has been shown that the unitary dynamics of quantum 
mechanics is the base of ensemble quantum computation [22, 23, 24]. The 
basis principle behind the ensemble quantum computation is that both a 
closed quantum system and its ensemble obey the same unitary dynamics 
of quantum mechanics if there is not any decoherence in both the quantum 
system and its ensemble. This basis principle allows one to use quantum en- 
sembles such as the spin ensembles without any quantum entanglement to do 
real quantum computation. Recently, the NMR multiple-quantum measure- 
ment techniques have been used to output quantum computational results 
which are the inphase multiple-quantum coherence spectra in a spin ensemble 
[24]. The NMR measurement for the inphase multiple-quantum coherences 
need not an exponential resource when the multiple-quantum coherences in 
the spin ensemble are created efficiently by any quantum circuit. Therefore, 
both the unitary dynamics of quantum mechanics and the multiple-quantum 
measurement techniques may form the base for the scalable ensemble quan- 
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turn computation [24]. In this paper both the unitary dynamics of quantum 
mechanics and the muhiple-quantum measurement techniques have been ex- 
ploited to solve the prime factorization problem on a spin ensemble without 
any quantum entanglement. The purpose for it is to study how the quan- 
tum dynamics plays an important role on the origin of power of quantum 
computation in the prime factorization on a spin ensemble. 



2. Eigenvalues and eigenvectors of modular exponential 

operation 

The integer factoring problem can be reduced to the order-finding prob- 
lem, while the latter is closely related to the unitary transformation of the 
modular exponential operation in the quantum factoring algorithms based 
on pure quantum states [5, 6]: 

U{y,r,N)\x) = \xymodN), x = 0, L ...,N- 1. (1) 
The transformation U {y, r, N) is a unitary transformation only when the 
numbers y and A'^ are coprime to each other. The explicit form of the uni- 
tary transformation is dependent only on the numbers y and A'^. This uni- 
tary transformation hides the period r of the modular exponential function 
f{m) = f{y,m,N) = mod that is a periodic function: /(m) = f{m+r). 
The period r need to be determined in the order-finding problem. The mod- 
ular exponential unitary transformation (1) can be efficiently implementable 
[1, 2, 5, 6]. Given a number y prime to the integer N the unitary transforma- 
tion U{y,r,N) can be determined explicitly from Eq.(l). Then the Hamil- 
tonian corresponding to the unitary operator U {y, r, N) can be expressed in 
form 

H{y,r,N)^i\nU{y,r,N). (2) 
Since the order of the unitary operator U{y,r,N) is r, i.e., U{y,r,NY — E 
(the unity operator) there are r different eigenvalues for the unitary operator: 
Afc = exp(— ■i27rA;/r), /c = 0, 1, r — 1. Then the Hamiltonian H{y, r, N) also 
have r different eigenvalues: = 27ik/r, k = 0,1, ...,r — 1. According to the 
Cayley-Hamilton theorem of linear algebra [25] the Hamiltonian H{y,r,N) 
of Eq.(2) can be expanded as 

Hiy,r,N)=f: a,U{y,r,N)K (3) 

fc=0 

Now suppose that the common eigenvectors of the unitary operator and its 
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Hamiltonian are denoted as {j^'jt)}, then their eigen-equations are respec- 
tively given by 

t/(y,r,7V)|*fe) =Afe|*fe) (4a) 

and 

H{y,r,N)\^k)=\k\'i^k). (46) 
By the operator equation (3) and the eigenvector one obtains from 
Eq.(4b) 

■r-1 

Afe =^ a; exp(— i27rA;Z/r). (5a) 

1=0 

Obviously, both the eigenvalues {A^ = 2'Kk/r} and the coefficients {ctfc} form 
a pair of Fourier transform, and according to the Fourier transform relation 
(5a) one can determine explicitly the coefficients {ak}, 

aik=i: ^exp(i27rM/r). (56) 
Then inserting Eq.(5b) into Eq.(3) the Hamiltonian H{y,r,N) is written as 

H{y, r, N) = E E ^ exp{i27Tkl/r)Uiy, r, N)K (6) 

k=Ol=0 

The Hamiltonian H{y, r, N) above is derived in detailed from its correspond- 
ing unitary operator U {y, r, N) of Eq.(l) partly due to that the manipulation 
for a Hamiltonian is usually more convenient than for a unitary operator in 
a complex spin ensemble [19, 20]. There is an important property for the 
Hamiltonian H{y,r,N) and unitary operator U{y,r,N) of the modular ex- 
ponential operation (1) according to the definitions of the unitary operator 
(1) and its Hamiltonian (2) that the unitary operator U{y,r, N)"^ can be 
written as 

U{y, r, N)"^ = Uiy"^, r, N) = exp[-imH{y, r, N)]. (7) 
The first equality in Eq.(7) shows that quantum circuit of the unitary oper- 
ation U{y, r, N)"^ can be efficiently constructed even when the integer m is a 
huge number, e.g., m = r, while the second equality shows that the integer 
m really acts as the discrete time variable in the dynamical process of the 
modular exponential operation. 

It is usually convenient to calculate time evolution of a spin ensemble un- 
der a spin Hamiltonian if the eigenvectors and their eigenvalues of the spin 
Hamiltonian are determined. The common eigenvectors of the unitary oper- 
ator U{y, r, N) and its Hamiltonian H{y, r, N) can be constructed explicitly 
using the eigenequations (4a) and (4b) . They should be a linear combination 
of the conventional computational base Since the period is r for the 
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modular exponential function /(m) = f{y,m,N) = y^modN, that is, there 
is the lowest integer r such that f{y,r,N) = 1, it follows from the unitary 
transformation (1) that for a given integer x one has 

U{y,r,N)\xy"' mod N) = \{xy^'^+^^'^°'^'-) modN), (8) 
where m — 0,1,. ..,r — 1. It is clear that the basis subset {\xy'' modN); 
A; = 0, l,...,r — 1} form a closed state subset S{x) under the unitary trans- 
formation U{y, r, N). For convenience, the dimension of the state subset S{x) 
is denoted as since it may depend on the integer x. Obviously, for x = 1 the 
dimension Tx of the subset S{x) equals r exactly. The dimension Tx is always 
smaller than or equal to the period r for any integer x : < x < N — 1, that 
is, the period r is the maximum dimension in the subsets S{x) for all possible 
x, 0<a;<A^ — 1 also because according to the definition (1) the unitary op- 
erator U {y, r, N) has the order r and is independent of any x, 0<a;<A^ — 1. 
In particular, Vx = 1 when x = and its subset >S'(0) = {|0)}. It follows from 
Eq.(8) that an arbitrary eigenstate |\E'<j(x)) of the unitary operator U {y, r, N) 
of the subset S{x) can be expressed as a hnear combination of the basis of 
the subset S{x), 

I'^six)) = f: c{rx,s,k)\xy'' mod N),s^ 0,1,.., Tx-l. (9) 

fc=0 

Inserting the eigenstate (9) into the eigenequation (4a) one obtains 
^ exp{— i2Tis/r)c{rx, s,k)\xy^ mod N) 

k=0 

r —1 

- J2 c{rx, s,k)\xy''+'^ mod N) ^0. (10) 

k=0 

If dimension r^, of the subset S{x) is r exactly then all r states {\xy'^ mod N) , 
m = 0,l,...,r— 1} are independent of each other. Then the recursive relations 

for the coefficients c(r, s, k) can be set up by the eigenequation (10) 

exp(— i27rs/r)c(r, s, k) — c{r, s, k — 1), k — 1,2, r — 1; 
and 

exp{—i27is/r)c{r, s, 0) = c(r, s,r — 1). 
Therefore, the coefficients c(r, s, k) are determined by 

c(r, s, k) = exp{i27isk/r)c{r, s, 0), k = 1, 2, r — 1. (11a) 
With the help of the orthonormal relations for the eigenstate |^s(a^)) and the 
basis \xy'^modN): ("^ s{x)\'$ s{x)) = 1 and (xy'^ modN\xy^ modiV) = 6kk' 
it is easy and straightforward to find the coefficient c(r, s, 0) = (here 
c(r, s, 0) is taken as a real). Then the other coefficients are given explicitly 
by Eq.(lla) once the coefficient c(r, s, 0) is known. Generally the modular 
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exponential function f{x, y, m, N) = xy'^ mod N {x ^ 0) may have the same 
period r as the function f{y,m,N) = y"^modN, but besides the period r 
the function f{x,y,m,N) may also have other periods different from r 
for some given x [14]. Then in the case of r^; < r the eigenstate subset 
{|\l's(a;))} of Eq.(9) of the unitary operator U{y,r,N) has both two period 
Tj. and r, and the period divides r because there must be the relations: 
\^q{x)) = |\['r(a;)) = |\E';r^(a;)), / is some integer. The eigenvalues of the 
unitary operator U{y,r,N) belonging to the eigenstates of Eq.(9) can be 
obtained using the eigenequation (4a), 

Uiy,r,NY^\-^s{x)) = A:^{x)\-^s{x)) = \-^s{x)), 
where the second equality is due to the period r-c, that is, Ixy'^modN) — 
Ixy''^'"'' modN). Therefore, the eigenvalues are given by 

As{x) = exp{-i27Ts/r^), s = 0,1, ...,r^ - 1. 
Again using the eigenequation (4a) and the eigenvalues A.s{x) one can set up 
the recursive relations for the coefficients c{rx-i s, k) of the eigenstate |\E's(a;)) 
similar to Eq.(lla) and hence the coefficients c{rx,s,k) are determined 

c(r^., s, k) = exp{i27isk/rx), k,s = 0, 1, - 1. (116) 
Then by using the coefficients of Eq.(lla) and Eq.(llb) one obtains from 
Eq.(9) the common eigenstates of the unitary operator U {y, r, N) and Hamil- 
tonian H{y,r,N) [A, 6, 7], 

\^s{.x)) = ^'Y. exp(?27rsA;/r,)|a;|/*^mod7V),s = 0,l,..,r^-l. (12a) 

fe=0 

Obviously, the eigenstate set {|^s(ic))} and the basis set {\xy^m.odN)} form 
a pair of Fourier transforms, and the inverse Fourier transform of Eq.(12a) 
generates the basis \xy^m.odN) as 

\xy^TiiodN) = ^ J2 exp{-i27rsk/rx)\^s{x)) (126) 
^''"^ fe=0 

The Fourier transforms of Eq.(12a) and (12b) are helpful for calculating in 
an analytical form the time evolution of a spin ensemble under the unitary 
operation U{y,r,N). This can be seen in next sections. 

Since every integer k in < k < N always can be expressed as A; = 
xy"^ mod N by choosing suitably the integers x and m, where the integer y 
is coprime to and smaller than the integer N, then the conventional com- 
putational basis can be expressed either as {|A:), k = 0,1, 2, — 1} or as 
{\xy"^ mod N) , m = 0,1, 2, r^, — 1; a; = 0, 1, N — 1}. One can classify 
the conventional computational basis {\k)} or {\xy™- mod N)} according to 
the transformation property of the unitary operator U {y, r, N). For example. 
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for a given integer x one can generate a basis subset S{x) — {\xy"^ modN), 
m = 0, 1, Tx- — 1}. The whole Hilbert state space with dimension N then is 
divided into t independent and orthogonal basis subsets S{x) with different 
integers x : xq < xi < ... < Xt-i. Clearly, = rxQ + rxj^ + ...+rx^_^. Therefore, 
the conventional computational basis set also can be expressed in the simpler 
form {\xiy''modN), k = 0,1, r-^j — 1; Z = 0, 1, t — 1} according to the 
transformation property of the unitary operator U{y,r, N) (1), which is also 
equivalent to the conventional computational basis set In particular, 

for xq = the subset S{xq) = {|0)} with = 1 and for xi = 1 the subset 
S{xi) is an r— dimensional subset. Actually, besides the subset S{xi) there 
may be also other r— dimensional subsets S{x) with x > 1. Suppose that 
there are d independent r— dimensional subsets including xi = 1 and x > 1 
in the A^— dimensional Hilbert space. There are rd computational basis that 
belong to the r— dimensional subsets among the N computational basis, and 
the rest N — rd basis are of those subsets with dimensions rx smaller than 
r. How many the computational basis belong to the r— dimensional subsets 
in all computational basis? This can be answered by the theorem (Parker 
and Plenio [14]): Given two prime numbers p and q, N = pq, r is defined 
as the period of the modular exponential function /(m) = f{y,m,N) = 
y"^modN for an arbitrary integer y, then there are at least {p — l){q — 1) 
positive integers x less than and coprime to the integer N such that the 
modular exponential function g^m) = f{x,y,m,N) — xy"^modN has the 
minimum period equal to r for < y < — 1. 

Actually, the unitary transformation (1) shows that the conventional com- 
putational basis set is also equivalent to the basis set {\xy'' mod N) , 
x = 0,l,...,A^ — l}for any given integer k. Then according to the theo- 
rem [14] number of the computational basis that satisfy xy^^^^ = xy^ mod A^ 
with r-r < r is (p + g — 1) at most in the whole basis set {\xy^m.odN), 
X = 0, 1, A^ — 1}. This also means that number of the computational ba- 
sis that belong to those basis sets S{x) with dimensions r^; < r is at most 
{p -\- q — 1) in the complete basis set {\xiy^ mod N) , k — 0,l,...,rxi — 1; 
I — 0,1, ...,t ~ 1}. Therefore, the number N — rd oi the computational basis 
of those subsets with dimensions r^, smaller than r is at most {p + q — 1), 
and the number rd of the computational basis of the r— dimensional subsets 
in the A?^— dimensional (A^ = pq) Hilbert space is at least pq — {p + q — 1). 

Finally, it is also important to know the orthonormal relations for the 
eigenstates \^s{,Xk)) and the conventional computational basis \xiy^modN) 
for conveniently calculating the time evolution of a spin ensemble under the 
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unitary operation U {y, r, N), 

{<l'kixi)\^k'{xi,)) = 6kk'5ii', (13a) 
{xiy'' mod N\x I' y''' mod N) = 5kk'5u'. (136) 

where k, k' = 0, 1, — 1 and 1,1' = 0,1, ...,t — 1. 



3. Time evolution process of modular exponential operation 

Generally the prime number N is not equal to some power of two. Sup- 
pose that the prime number N = pq satisfies 2"^"^ < N < 2^. For simplifying 
calculation of the time evolution of a spin ensemble during the modular expo- 
nential operation the NMR quantum computer could be chosen conveniently 
as a heteronuclear spin ensemble Iil2---IniSiS2---Sn (denoted briefly as ImSn) 
that consists of spin- 1/2 I nuclei and n spin- 1/2 S nuclei and particularly 
Tii = 1,2,.... This is just like the pure-state Shor's factoring algorithm us- 
ing two memories [1, 2]. The conditional modular exponential operation 
Ui^sSv-i '^1 ^) s-PPlyi^g to the spin system In^Sn is built up with the unitary 
transformation of Eq.(l), 

\a){U{y,r,N)\x)), x = 0, 1, TV - 1; a = l 
Ui,sSy,r,N)\a)\x)^{ \a)\x), x = 0,1, ...,N - l; a = (14) 

\a)\x), X = N,N - I; a = 0,1 

where the quantum states \a) and belong to the kih spin / and all n spins 
Sn of the spin system ImSn, respectively. The modular exponential operation 
U {y, r, N) is applied only to those quantum state |a;) of the spin subsystem Sn 
with X < N only if the kth 7— spin quantum state \a) = |1). It can turn out 
that the Hamiltonian corresponding to the unitary operator Ui^sSu-if-i^) 
(14) may be expressed as 

Hi^sSy.r,N) = E^^...^Ej,_^ 

®{\E,-h,) (g) (g) ...Ki (g) HsSy. r, N) (15) 
where /fe|0) = ||0) {h = 1) and Ikz\l) = — and the Hamiltonian 
Hs^ {y, r, N) applied only to the subsystem Sn is defined as 

Hs,, {y, r, N) = H{y, r, N) Z^^^ (16) 
where Zl_]^ is the (L — N) x (L — iV) — dimensional zero operator. It is easy to 
prove according to Eq.(7) that the conditional modular exponential operation 
C//^5^(y, r, A^)"^ can be expressed as 

Uj,sAy,r,Nr = Uj,sM\r,N) = exp[-M45„(y, r, iV)]. (17) 
The unitary operator Ui^s„{y,i^, N)"^ can be constructed efficiently for any 
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integer m because the unitary operator Uii^s„iy"^,f,N) can be constructed 
efficiently. It is known from Eq.(17) that the power m in the unitary operator 
Uj^SniU) ^1 N)™- really acts as the discrete time variable in the dynamical pro- 
cess of the conditional modular exponential operation. Note that any pair 
of the Hamiltonians Hi^sSv-,^-,^) (15) with different I-spin operators Ikz 
{k — 1,2, ...jTii) commute each other. A more general conditional modular 

exponential operation can be constructed by 

m 

Uir.,sAy,r,N)=U Ui,sMr,N), (18) 
k=i 

and it is easy to prove that the relation (17) also is met for the general condi- 
tional modular exponential operation (18). The Hamiltonian of the unitary 

operator Ui^.s„{y,r, N) therefore is written as 

m 

Hr^^sAy,r,N) =E Hr^s„{y,r, N). (19) 
k=i 

Obviously, the conditional unitary operation CZ/^.^^ {y, r, N) is independent 
of any quantum state |a)|a;) (as the initial input state) of the spin system 
ImSn- This suggests that the unitary operator can be applied not only to 
any pure quantum states of the spin system /„^S'„ but also directly to any 
mixed states of the spin ensemble ImSn of the spin system [22, 23]. This is 
just the essence of the basic principle that both a closed quantum system and 
its ensemble obey the same unitary dynamics of quantum mechanics if there 
is not any decoherence effect in both the quantum system and its ensemble 
[22, 23, 24]. This principle forms the base of the current factoring algorithm 
and the real implementation of the algorithm on a spin ensemble. The uni- 
tary operator Uj^^s^iy, f-i N) hides the period r to be determined. In order to 
find the period it ffi^st needs to transfer the information of the period r of the 
unitary operator into quantum states of a quantum system or the density op- 
erator of its quantum ensemble because both the quantum states and density 
operators can be measured conveniently in practice. According to the ba- 
sic principle the initial input state of the unitary operator C//^5„(z/, r, N) can 
take either any pure quantum state of a closed quantum system or any mixed 
state, i.e., density operator of its ensemble, but in an NMR spin ensemble it 
is most convenient to take the initial density operator, i.e., the input state of 
the current factoring algorithm, as the thermal equilibrium state of the spin 
ensemble. For the spin ensemble ImSn in a high magnetic field the thermal 
equilibrium state can be written as, in high temperature approximation 

ni n 

p,^ = aE+ ^^kh.+ E ^skSkz (20) 

k=l k=l 
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where E is the unity operator and the operators Ikz and Skz are the longitu- 
dinal magnetization operators of the kth spins / and S, respectively. Then 
a nonselective 90° excitation pulse i?j(90°) = exp(— i7r/j^/2) applied to all 
the spins / and a 90° nonselective pulse Rs{90'^y) = exp[=pi(7r/2)5'y] with 
two-step phase cycling if — +y, —y applied to all the spins S convert the 

thermal equilibrium state (20) into the single-quantum density operator, 

m 

P(0) = (E eikIkx)<S)Et<S)EI<S) -<S)K (21) 
fe=i 

where the unity operator term aE is neglected without losing generality and 
E^ is the 4— dimensional unity operator of the kth spin S. The two-step 
phase cycling (p — +y, —y [19, 20] cancels the contribution of the thermal 

n 

equilibrium magnetization p^g^ = ^ SskSkz of the subensemble Sn of the spin 

k=l 

ensemble ImSn to the output NMR signal, leaving only the thermal equihb- 

ni 

rium magnetization p^^^ =J2 ^ikhz of the subensemble In, having a net 

fc=i 

contribution to the output NMR signal. Now the information of the period 
r in the unitary operator [//^_5^(y, r, iV) need to be loaded on the density 
operator of the spin ensemble ImSn- This can be achieved by applying the 
unitary operator on the initial density operator p(0). An analytical calcula- 
tion is important for the time evolution of the spin ensemble with the initial 
density operator p(0) under the conditional unitary operation Ui^,s^{y, r, N). 
It can be performed conveniently by first expressing the initial density op- 
erator (21) in terms of the common eigenvectors of the unitary operator 
U{y,r. N) (1) and its Hamiltonian H{y,r, N) (2). By using the conventional 
computational basis set {\xiy^ modN), k = 0,1, — 1; / = 0, 1, t — l) 
the initial density operator (21) is rewritten as 

rii t—V'^i^^ 

P(0) = (E ^ikhx)^ E E \xiy^m.odN){xiy^m.odN\ 

k=l 1=0 k=0 

ni L—1 

+ {E^^kh■■c)^ E \ks){ksl {L = 2'\\L<N <L). (22) 

k=l ks=N 

With the help of the inverse Fourier transform of Eq.(12b) and the orthonor- 
mal relation (13b) the density operator p(0) is further expressed as 

p(o) = (E^ife4.)(8)E E |*.(xfe))(*.(xfe)| 

fc=l k=0 s=0 

ni L—1 

+(E £ife4a.)(8) E (23) 

k=l ks=N 
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On the other hand, with the help of the definition (14) of the conditional 
unitary operation [//^^^(i/, r, A^) and its the Hamiltonian (15) as well as the 
eigenequation (4b) of the Hamiltonian H{y,r,N) it is now easy to calcu- 
late the time evolution of the spin ensemble when applying the conditional 
unitary operation Ui^,s„{y,r, N)™- on the density operator p(0) of Eq.(23), 

pM = Ui^^sAy,r,Nrp{o)ui^sAy,r,Nr 

^-Eily^Y. E sm[2Tims/r^^]\^!s{xi)){^!s{xi)\ 

l=Q s=0 

t-Vxi-l 

+£iIx^E E cos[27rms/r^^]\'^s{xi)){'^s{xi)\ 

1=0 s=0 

+eih® Z \ks){h\, (24) 

ks=N 

ni 

where all rij spins / have the same spin polarization factor £j, that is, £j/^ = ^ 

k=l 

^iklkn ifJ' = x,y,z). One can obtain the antisymmetric y— component of the 
density operator (24), 

Py{m) = -eiIy<S)E E sin[27rms/r,J|*,(xO)(*s(x;)|, (25) 

1=0 .8 = 

by doing another experiment: p{—m) = t/^^ g^(|/, r, A^)™p(0)t//^ ,s„(|/, r, A^)™ 
and then coadding coherently the final output NMR signals of the two density 
operators p{m) and [— p(— m)], that is, Pyim) = ^(p(m) — p{—m)). The 
inverse unitary operation Uj'^ SniV- ^' can be also implemented efficiently 
just like the unitary operation f//„,5„(y, r, N)"\ as can be seen in next section. 
The density operator Py{m) of Eq.(25) is antisymmetric for any integer m : 
py{kr /2 — m) = —py{kr/2 + m), k — 0,1,2,..., because the period r can 
be divided by any dimensions r^.^. By inserting the eigenstates \^ s{xi)) of 
Eq.(12a) into Eqs.(24) and (25) one can express the density operators p{m) 
and Py{m) in terms of the conventional computational basis, 

p{m) = -Sily^'Y. E {^;}{-\xiy^TiiodN){xiy^+^uiodN\ 

1=0 k=0 

E E {\){\xiy'' mod N){xiy^+''' mod N\ 

1=0 k=0 

+ \xiy^+'^modN){xiy''modN\}+e,I,^ ^ \K){K\ (26) 

ks=N 

and 
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Py{m) = -Eily^ Yl E {^){-\xiy^ mod N){xiy''+"' mod N\ 

1=0 k=0 

+ \xiy''+"' mod N){xiy'' mod N\}. (27) 
Note that the conventional computational basis {{xiy'^ mod A^), / = 0, 1, 
t — l;k = 0,l,...,ra;, — 1} of the subensemble Sn are orthogonal to each other, 
as shown in Eq.(13b). Then the operator \xiy'^ mod N){xiy'''^'^ mod N\ is a 
diagonal operator only when the equality xiy'' mod N = xiy^^'^ mod N holds. 
As shown in the previous section, the equahty holds only when m = k'r^i 
{k' = 0,1,...,). Therefore, the operator |a;/y^modiV)(a;;y*^"'"'"modiV| is an 
off-diagonal operator when m ^ k'vxi- It is known that the diagonal el- 
ements of a density operator are the noncoherence components which are 
known as the conventional longitudinal magnetization and spin order com- 
ponents in NMR spectroscopy, while the off-diagonal elements represent the 
coherent components of the density operator, and it is also well known in 
NMR spectroscopy that the coherent components of a density operator are 
the conventional multiple-quantum coherences including single-quantum co- 
herence [19, 20]. It is shown below that the density operator Py{m) of Eq.(27) 
is a pure multiple-quantum coherence operator of the subensemble Sn- Here 
assume that the period r is an even integer. First, the density operator py{m) 
is clearly a pure multiple- quantum coherences when m ^ k'rx^- Next, it need 
to be shown that the integers m — k'r^i are the zero points of the density 
operator Py{m). The zero points of the density operator Py{m) are defined 
as those integers m satisfying Py{m) = 0. Obviously, m = k'r are the zero 
points of the density operator Py{m) because the period r can be divided by 
any dimensions r^i and hence m = k'r = k^iVxi {kxi is an integer) which lead 
to the identity: 

rx,-l 

^ Wxiy'^^"' mod N){xiy^ mod N\ 

k=0 

-\xiy'' mod N){xiy^+''' mod N\] = 0, (/ = 0, l,...,t-l), (28) 
and thus, the density operator Pyirn) = 0. Moveover, it is easy to prove 
that every integer m = {2k' + l)r/2 satisfies the formula (28) and hence 
m — {2k' + l)r/2 are also the zero points of the density operator Py{m). 
Therefore, all possible zero points of the density operator Py{m) are given 
by m = k'r/2 (A:' = 0, 1, 2, ...). Then for any integer m the density operator 
Py{m) is either an off-diagonal operator or equal to zero. This indicates that 
the density operator Py{m) is a pure multiple-quantum coherence operator 
of the subensemble Sn and does not contain any longitudinal magnetization 
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and spin order (LOMSO) operators of the spin subensemble. 

The density operator p(m) of Eq.(26) is more comphcated and consists of 
both the LOMSO operators and multiple-quantum coherence operators of 
the subensemble Sn- The first term on the right-hand side of Eq.(26) is re- 
ally the pure multiple-quantum coherence operator Py{m), while the last term 
really consists of pure LOMSO operators. These LOMSO operators are re- 
ally invariants under the unitary transformation in Eq. (24) with the unitary 
operator ?7/,j.5',j(y, r, A^)™". The second term contains the invariant diagonal 
element Poo — which keeps unchanged under the unitary transforma- 

tion. When the integers m = k'r the density operator of Eq.(26) equals the 
initial density operator: p(m) = £^7^ = p(0), indicating that at the points 
m = k'r the initial density operator p(0) is not transferred into any multiple- 
quantum coherences of the subensemble Sn under the unitary transformation. 
Actually, if m = k'r the unitary operator Ui^,s„{y,r, N)"^ is the unity op- 
erator, that is, Ui^,s^{y,r, N)"^ = E, and any initial density operator keeps 
unchanged by the unity operation. One sees that at m = k'r there is not 
any multiple-quantum coherence of the subensemble Sn in the density oper- 
ator p{m) and also in the density operator Py{m). Therefore, these integers 
m = k'r ik' = 0,1,...,) are really the zero points of the density operators p{m) 
and Py{m) in the sense that there is not any multiple- quantum coherence of 
the subensemble Sn in these density operators, although p{m) — p(0) 7^ at 
m — kr. When m = kr,j.^ the whole conventional computational basis subset 
S{xi) = {\xiy'' mod N) , k = 0,1, ...,r^; — 1} of the subensemble Sn keeps un- 
changed under the unitary operation Ui^,s„{y,r, N)"\ showing that the rxi 

operator terms (J^ ^ikhx) ^ ^ \xiy'' mod N){xiy'' mod N\ of the initial 

k=l k=0 

density operator p(0) of Eq.(22) can not be transferred into the multiple- 
quantum coherences. Based on these facts one can calculate the conversion 
efficiency of the initial density operator p(0) into the multiple-quantum co- 
herences under the unitary transformation. 

It follows from Eq.(21) that the initial density operator p(0) consists of 

ni 

L operator terms p^^ = Siklkx) ^ \ks){ks\, kg = 0, 1, ...,L — 1. That the 

k=l 

density operator p(0) can not be transferred completely into the multiple- 
quantum coherences of the subensemble Sn is because the L — N opera- 

L— 1 Til 

tor terms Pks ^ V2) and = (S ^ife4a;)|0)(0| of the den- 

ks=N k=l 
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sity operator p(0) keeps unchanged under the conditional unitary operation 
f//^ r, A^)™" with any integer m. But the total contribution of these in- 
variant operator terms to the density operator p(0) is clearly less than 50%. 
Therefore, when m ^ krxi {I = 0, 1, t — 1) the initial density operator p(0) 
is transferred into the multiple-quantum coherences in a high efficiency of 
50% at least, and this efficiency is independent of the qubit number n of the 
spin ensemble {IrnSn}- For the case m = kr^^ (each r^^ < r) but m ^ k'r the 

ni 

operator terms (^2 ^ikhx)® Yl, \xiy^ Ta.odN){xiy^ m.odN\ also keep 
unchanged under the conditional unitary operation in addition to the L — N 

L-l 

operator terms J2 Pk^ Poo^ t)ut the contribution from all these invari- 

ks=N 

ant operator terms to the initial density operator p(0) is at most {L — rd)/L, 
as shown in the previous section. Then the conversion efficiency of the ini- 
tial density operator p(0) into the multiple-quantum coherences under the 
conditional unitary operation will be rd/L — {N — p — q + 1)/L at least 
[14], as can be seen in the previous section. Note that < N < L. The 
efficiency rd/L is generally not less than 50% for a sufficient large number 
N = pq, indicating that for an arbitrary m ^ kr the initial density operator 
p(0) is efficiently transferred into the multiple-quantum coherences of the 
subensemble Sn in the density operator p{m) with an efficiency generally not 
less than 50% under the conditional unitary operation Ui^,s„{yTf,N)"^ on 
the spin ensemble ImSn with a large qubit number n. This high conversion 
efficiency will directly result in that the multiple-quantum coherences of the 
subensemble Sn in the density operator p(m) may be detected efficiently, as 
can be seen in next section. 



4. NMR measurement of multiple-quantum coherences 

Generally all the multiple-quantum coherences can not be observed di- 
rectly in the NMR measurement except the single quantum coherence. All 
non-first-order multiple-quantum coherences usually may be detected indi- 
rectly through the direct NMR measurement of the single-quantum coher- 
ence. To detect multiple-quantum coherences one first needs to convert them 
into single-quantum coherence by making a unitary transformation composed 
of a sequence of RF pulses and interaction intervals in the spin ensemble. In 
general, the density operator of the spin ensemble, for example, p{m), can 
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be expanded as a quantum coherence order series [19], 

n 

p{m) = E (7p(m), (29) 

p=—n 

where the operator ap{m) is the p— order quantum coherence operator and 
(Tp{m) = a-p{m)~^ because the density operator is an Hermitian operator. 
Since the present discussion is focused on the multiple-quantum coherences 
of the spin subensemble Sn of the spin ensemble In.Sn the order index p 
of the expansion series (29) in this special case is referred to the p— order 
quantum coherence of the spin subensemble Sn- The maximum quantum 
order of the subensemble Sn is n. This means that the density operator 
p{m) (29) are generally composed of the multiple-quantum coherences with 
at most 2n + 1 different quantum orders range from —n to n. The density 
operator p(m) can be detected through the spins / or S* by applying a unitary 
transformation to convert it into single quantum coherence. Here consider 
the measurement method through detecting the single quantum coherence of 
the spin / instead of the spin S. Before the multiple-quantum coherences are 
transferred into the single-quantum coherence of the spin / they are labelled 
with their own precession frequencies in order to distinguish different order 
quantum coherences, and in order to observe effectively the multiple-quantum 
spectral peaks of the multiple-quantum coherences one had better label all 
the same order quantum coherence with a single precession frequency. Then 
the spin Hamiltonian used to label the multiple-quantum coherences of the 
subensemble Sn may be chosen as 

n 

Hs ^Yl ^SkSkz = ^^sSz- (30) 

k=l 

This labelling Hamiltonian is independent of the spins /„.. Under this Hamil- 
tonian the multiple-quantum spectrum of the spins Sn has at most 2n + 1 
different order multiple-quantum peaks with their own precession frequen- 
cies, and the frequency for all p— order quantum coherence is simply equal to 
P^s ip — —n, —n -|-l,...,n — l,n). The frequency labelling for the multiple- 
quantum coherences may be achieved by the time evolution process of the 
spin ensemble (/„.5'„) starting at the density operator p{m) under the Hamil- 
tonian (30), 

n 

p{7n,ti) = exp{-iHsti)p{m)exp{iHsti) = ap{m) exp{-ipujsti). (31) 

p=—n 

Then the density operator (31) is converted into single-quantum coherence 
under the specific unitary transformation V/^_5^(i/, r, iV) which may be con- 
structed through the modular exponential operation U{y,r,N), 
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p/m,ti) = Vi^^^sSy.r,N)p{mM)Vu^sMr,N)+. (32) 
Then the observable NMR signal for the density operator pf{m,ti) is given 
by 

Sf{m,ti) = Tr{Fpf{m,ti)} 

= Tr{Vj^My, r, N)+FVj^^sM r, N)pf{m, h)}, (33) 
where the observable single-quantum operator F — I^, which is proportional 
to the initial density operator p(0) (21). Because all the same order quantum 
coherence has the same precession frequency the total amplitude for a given 
order quantum NMR signal is the coherent sum of the amplitudes of all the 
same order quantum coherence. The unitary operator Vi„,s„{y,r, N) should 
be constructed suitably in order that the coherent sum of the amplitudes 
is constructive and the inphase multiple-quantum spectrum is generated for 
all the same order quantum coherence, otherwise the total amplitude may 
severely attenuate due to the destructive coherent sum. One of the best 
ways to build up the unitary operator is simple to make the unitary operator 
satisfy [26], 

Vr^^sAy,r,N) = Ul^sAy,r,Nr. (34) 
This is a direct requirement of the famous time- reversal symmetry [27]. In 
high-resolution NMR spectroscopy the time-reversal symmetric unitary oper- 
ator usually may not be easily created from the scalar J-coupling interactions 
of a complex coupled spin system in a liquid. An alternative method to gen- 
erate the inphase multiple-quantum spectrum for all the same order quantum 
coherence may be that a series of experiments are performed starting from the 
same density operator, e.g., p(m, ti) of Eq.(31), and using different unitary 
operators V}^^g^{y,r,N) so that the following relation is satisfied: 

E ^) = ^rV(m)+ = E £rVp(m)+, (35) 

k ' ' p=—n 

then the inphase multiple-quantum spectra could be generated as well by 
adding coherently these experimental NMR signals, although the Hamilto- 
nian of every unitary operator V^^\^{y^r^N) may not always satisfy the 
time-reversal symmetry [27]. The inphase multiple-quantum spectrum of the 
density operator Py{m) of Eq.(25) may be generated by this method later. 

According to number theory [28] Euclid's algorithm may be used to find 
efficiently the multiplicative inverses in modular arithmetic. The inverse of 
the integer y modular A^, i.e., y~^modN can be efficiently calculated by 
Euclid's algorithm by taking 0((log2 A^)^) steps [28]. Then according to the 
definition (1) of the unitary operator U {y, r, N) the inverse unitary operator 
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U(y, r, N)'^ then can be constructed by 

U{y,r,N)+ = U{y-^modN,r,N). (36) 
Then the unitary operator U{y,r,N)^ can be efficiently constructed just as 
the unitary operator U {y, r, N). Consequently the inverse conditional unitary 
operation Ui^ sAy^ ^' '^^^ ^® constructed efficiently by the definition 
(14) of the conditional unitary operation Ui^,s^{y,r,N). This directly leads 
to that the unitary operator U^^ g^{y,r, N)'^ can be efficiently constructed 
for any integer m because U^^ g^iy, r, N)"^ = g^{y^, r, N). 

As an example, below the NMR multiple-quantum signal Sf{m,ti) of 
Eq.(33) is calculated in detailed using the unitary operator V/^_5^(t/, r, A^) of 
Eq.(34). There are the orthogonal relations between different order multiple- 
quantum coherence operators, 

Tr{ap(m)+(jp'(m)} = Tr{\ap{m)\'^}Spp'. (37) 
Inserting Eq.(31) and (34) into Eq.(33) and using the orthogonal relations 
(37) one can further write the NMR signal (33) as 

n 

Sf{m,h)^Y: e;'Tr{\ap{m)\^}exp{-tpu;sh). (38) 

p=—n 

The multiple-quantum spectrum Sf(m,U! — poos) can be obtained by fast 
Fourier transforming (ti) the time-domain NMR signal Sfirri^ti) of Eq.(38). 
Then the p— order quantum peak in the multiple-quantum spectrum has an 
intensity I{p,m) = e^^Tr{\ap{m)\'^} and the total intensity for the multiple- 
quantum spectrum is given by 

n 

/M= E e;'Tr{\apim)\'}. (39) 

p=—n 

The expansion (29) and the orthogonal relations (37) show that the total in- 
tensity I{m) also can be expressed as 

I{m) = e-^Tr{p{m)+p{m)}. (40) 
The formula (39) and (40) are really the direct result of equation. (35), a gen- 
eral method to obtain inphase multiple-quantum spectra. The total spec- 
tral power J(t) of the density operator p{t) = U{t)p{0)U{t)^ including the 
contributions from both the multiple-quantum coherence and the LOMSO 
operator components actually keeps unchanged when an arbitrary unitary 
operation U{t) is applied to the initial density operator p(0), 

m = e-'Tr{\pim = e-^Tr{\piO)n. (41) 
This property could be hclpftil for conveniently manipulating the modular 
exponential operation. Inserting all the multiple-quantum coherence compo- 
nents of the density operator p{m) of Eq.(26) into (40) the total intensity 
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/(m) can be calculated explicitly using the orthogonal relations (13a) and 
Tr{I^.y} — |ni2"' (the trace is only in the subensemble /„.), 

I{m) = ln,2"'e,{N- EE r,J{m,jr,J}, (42) 

j=Ok=0 

where the delta function S{x,y) = S^y The second term of Eq.(42) accounts 
for those invariant diagonal operators {\xiy'^ mod N){xiy'' mod N\} under the 
conditional unitary operation in Eq.(24) which have not a net contribution to 
the multiple-quantum spectrum. The total intensity /(m) satisfies, as shown 
in previous section, 

/(m) = if m = 0, r, 2r, 

I{m) > |ni2"'£j[iV - (p + g - 1)] if m 7^ 0, r, 2r, .... 
If one detects only the multiple-quantum coherences of the antisymmetric 
density operator Py{m) of Eq.(25) then the total intensity Iy{m) of the 
multiple-quantum spectrum of the density operator Py{m) can be derived 
as 

Iy{m) = er'Tr{py{m)+Py{m)}. (43) 
Since the density operator py{m) is a pure multiple- quantum coherence op- 
erator one can calculate the total intensity Iy{m) directly by inserting the 
density operator Py{m) of Eq.(25) into Eq.(43), 

Iy{m)^e;'Tr{\eiIylS)E E sm[2nms/r,,]\^,{xk)){^s{xkW} 

k=0 s=0 

= |n,2-'£,[iV- EE r.A'^m,jr,^)], (44) 

j=Ok=0 

where the orthogonal relations (13a) of the eigenvectors l^^si^k)) a-nd the re- 
lations below has been used, 

E exp(±i47rms/r^J = r^J{2m,jr^^),j = 0, 1,2, ... (45) 

s=0 

As can be seen in previous section, each r^^. can divide the period r and the 
total number of the computational basis of those basis sets S{xk) with the 

periods r^j. < r is {p + q — 1) at most including the state |0) [14], that is, 

t-i 

for any given integer m, {p + q — 1) >E E ^a;fc<^(2m, jr-E^.) with the sum for 

j=Ok=0 

index k running only over all r^;^ < r, and most of the computational basis 
(>A'^ — (p-|-g — 1)) belong to those subsets with the period r. Therefore, 
the total intensity Iy{m) satisfies, 
Iy{m) = if m = 0, r/2, r, 3r/2, 

lylm) > \ni2'''ei[N - (p + g - 1)] if m 7^ 0, r/2, r, 3r/2, .... 
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Here suppose that the period r is an even integer, otherwise r /2,2)r /2, .... are 
not integers and hence not zero points and the density operator Py{m) has 
only zero points with integer m = 0, r, 2r, .... 

If the initial thermal equilibrium density operator of Eq.(20) now is trans- 
ferred completely into the observable single quantum coherence, then the 

n 

generated NMR spectrum will have a total intensity Iq = £jTr{/f}+ ^ 

k=l 

^skTr{Sl^} = ^2^'^~^^*\ni€i + nes). The conversion efficiency from the ther- 
mal equilibrium state of Eq.(20) into the multiple-quantum coherences of the 
density operators p{m) and Py{m) under the conditional unitary operation 
in Eq.(24) are given by the ratios Iy{m)/ Iq and Iy{m)/ Iq, respectively, which 
satisfy, 

I{rn)/ Iq = if m = 0, r, 2r, 

/(m)//o > riiSilN - (p + q - l)]/[2^{niei + ne,)] if m 7^ 0, r, 2r, .... 
and 

Iy{m)/lQ = if m = 0,r/2,r, 3r/2,..., 

Iy{m)/lQ > luiSilN - ip + q-l)]/[2^in,e, + nes)] ifm^0,r/2,r, 3r/2,.... 
In fact, the multiple-quantum coherences of p(m) (26) and Py{m) (27) are 

generated only from the thermal equilibrium state p^^^ = ^ £iklkz of the spins 

k=l 

7„. (which is completely transferred into p(0) (22)), while the thermal equilib- 
rium state magnetization (p^eg) ^^e spins Sn of the initial thermal equilib- 
rium state (20) has not a net contribution to the multiple-quantum coher- 
ences. The total intensity for the multiple-quantum spectrum generated from 
the initial density operator p(0) (22) should be 7^0 = ^iTr{Il] = jniSi2^''+"'\ 
Therefore, the conversion efficiency for the initial density operator (22) into 
the multiple-quantum coherences Py{rn) should be measured more exactly 
by the ratio Iy{m)/IiQ instead of the ratio Iy{m)/lQ. The ratio Iy{m)/IiQ 
satisfies, 

Iy{m)/IiQ = 0, m = 0, r/2, r, 3r/2, 

Iy{m)/IiQ >^[N-{p + q- l)]/2", m ^ 0, r/2, r, 3r/2, .... 
It can be seen that the efficiency is almost independent of the qubit number 
n of the spin ensemble {In.Sn) for a large prime integer N = pq (2"~^ < 
N < 2"'). Therefore, for a large prime integer and m 7^ 0,r/2,r, 3r/2, ... 
the initial thermal equilibrium density operator p^^^ is efficiently transferred 
into the multiple-quantum coherences Py{m) with an efficiency more than 
1/4. Since there arc 2n + 1 peaks in the multiple-quantum spectrum of Py{m) 
the intensity for each peak of the 2n+l peaks, on average, is approximately 
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inversely proportional to the qubit number n, indicating that intensities for 
some of the 2n + 1 multiple-quantum peaks do not reduce exponentially as 
the qubit number n. 

The density operator Py{m) of Eq.(27) is a pure multiple-quantum co- 
herences of the spin subensemble Sm but p{m) contains both the multiple- 
quantum coherence and the LOMSO operators. The LOMSO compo- 
nents may hamper the detection of zero-quantum coherence in the multiple- 
quantum spectra since both have zero frequency in the multiple-quantum 
spectra under the frequency labelling Hamiltonian (30). A better method 
may be using only the density operator Py{m) and its intensity Iy{m) to 
solve the factoring problem. According to the scheme (35) to create inphase 
multiple-quantum spectra the pulse sequence to create the output NMR sig- 
nal with the intensity lyim) (43) consists of the two experiments: 
Pf{m,h) = UlsAy"'.r,N)expi-tHsh)pyim)exp{iHsh)Uj^^sAy''\r,N), 

xpyim) eMiHsti)Ri{lSOl)+Ul^sAy'^^ r, N). 
By adding coherently the output NMR signals of the two experiments one 
will obtain the desired NMR signal with the intensity Iy{m) (43). Since the 
density operator Py{m) is obtained from the two experiments, as shown in 
previous section, the complete pulse sequence to create the intensity Iy{m) 
consists of four experiments. 

A more general initial density operator p(0) is suggested below for the 
modular exponential operation sequence above. For convenient treatment, 
the initial density operator p(0) still has a general LOMSO operator P5^(0) 
of the subensemble Sn-, that is, 

p(0) = aE+ E P/„,(0),(8)P5„(0),. (46) 

The two components pi^ (0) and P5'„(0) of the density operator p(0) belong 
to the two subensembles /„. and Sn of the spin ensemble ImSn, respectively. 
Particularly, in previous factoring sequence the density operator component 
Ps„{Q) is the unity operator E, as can be seen in the initial density operator 
of Eq.(21). The initial density operator can be prepared properly from the 
thermal equilibrium state (20) of the spin ensemble /„-S'„ by a suitable pulse 
sequence. Ignoring the unity operator term aE the initial density operator 
p(0) (46) is generally written as 

ni t-l^xi-l 

P(0) =S I] Pj{xiy^ mod N)\xiy'' mod N){xiy'' mod N\ 

j=l 1=0 k=0 
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ni L—1 

+ E e,Ij.<S) E P,(^)l^)(fc|, (L = 2ML < TV < L). (47) 

j=l k=N 

where the coefficient Pj{xiy'^ mod N) is the diagonal element with the index 
{xiy'^ mod N) of the initial density operator component pg^{0)j. There are 
the unitary transformations according the definition (14) of the conditional 
modular exponential operation: 

Ui,s,. {y, r, N)^Iy, (g) \xiy^ mod N) {x^y^ mod Np+g^ {y, r, N)^ 
= |(|l)(0|)fc(g) |.Ti?/+™modiV)(,T/?/modA^| 

+\{\Q){l\)k®W/^odN){xiy^+"'modNliij = k; (48a) 

= ( 1 1 ) ( 1 1 ) fc (g) Jjo; (g) I mod AT) mod TV I 

+(|0) (0|)fe (g) Ijx (g) \xiy'' mod N) {xiy'' mod , if j ^ k. (486) 
The unitary transformations (48a) and (48b) show that the diagonal oper- 
ator mod A^)(x;|/*^modiV| of the subensemble Sn can not be converted 
into multiple-quantum coherences but into other LOMSO operators by the 
conditional modular exponential operation Uii^g^{y,r, N)"^ {k ^ j) which 
Hamiltonian (15) does not contain the operator Ij^. The unitary trans- 
formations can be used further to calculate the time evolution of the spin 
ensemble with the initial density operator p(0) (47) under the conditional 
modular exponential operation (18), 

pM = Uj^^sAy,r,Ny^p{0)Ulg^{y,r,Nr 

= EE E ejp^{xiy''modN)Ui.s^{y,r,Nr 

j=U=0 k=0 

x{(^ii)i (g) ... <S)iEn)j-i (g) I,. <S)iEu)j+i (g) ... (g)(^ii)n, 

(g) modN){xiy''+"<''^-^^ modN\ 

+ [(Eoo)l (g)(^ll)2 (g ... (g(^ll),-l (g Ij. (g(^ll),+l (g ... (g(^ll)n, 
+ {Eu)l (g(^00)2 (g(^ll)3 (g) ... (g)(^ll),-l (g) Ij. (g)(^ll),+l (g - 

(S){EnU + ■■■ + (^ii)i - <S){En)j-i (g Ijx (g(^ii)j+i (g - 
(g(£;ii)„._i (g(£^oo)nJ (g mod A^)(x,y'=+"*("*-2) mod 

+ [(^00)l '(g)(^00)2 (g)(^ll)3 (g) ... (g)(^ll),-l (g (g)(^ll),+l (g) - 
(g)(En)n, + ... + (^ll)l (g) ... (g)(^ll)i-l (g Ijx (g(^ll),+l (g) ... (g)(^ll)n,-2 

0(^oo)n,-i (g(^oo)nj (g mod A^) mod A^] 

+ ... + (Eoo)l (g)(^00)2 (g) ... (g)(^00)i-l (g) /jx 0(£^oo)i+l (g ... (g)(^00)ni 

(g \xiy'' mod AT) mod N\}U^.s^ [y, r, N)"" 

ni L—1 

+ E^./..(g) E P,(fc)|A:)(A:|, (49) 
i=i fc=Af 

where the operator E^j — and the 2 x 2— dimensional unity operator 
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of the kth spin / is expressed as El = (i?oo)fc + (-E'ii)ik- The density operator 
p{m) (49) can be simphfied by using the shift-invariance identity: 

E E Pjixiy'' mod N)\xiy''+'i"' mod N){xiy''+i"' mod N\ 

1=0 k=0 
t-irxi-l 

= E E Pjixiy^''-'^"''^'^°'^''^i mod N)\xiy'' mod N){xiy'' mod N\, (50) 
;=o k=o 

where < (A; — qm) modrj;; < for any integers q and m. For example, 
one of those operator terms in the density operator p{m) (49) is calculated 
in detailed below, 

p«(m) = Ui^sAy:r,Nr{Enh0...0{En)j-i0Ij.0{En)j+i 

<S> - <S>{Eii)n, <S) E E ejp^ixiy'^ mod N) 

1=0 k=0 

X \xiy''+i'^ mod N) {xiy''+i'" mod N\Uf.s^ {y, r, N)"^ 

= |(^ii)i (8) ... (8)(£^ii),-i (8)(|i)(o|)/(8)(^ii),+i (8) ... (8)(^ii)n. 

(8) E E eiPj(a:ii/[*^-«™]'"°'i''-i modiV)|x/y^'+™modiV)(,Ti|/^modAr| 

+i(^n)"i°(8) ... (8)(^ii).-i (8)(|o)(i|),- (8)(^ii),+i (8) ... (8)(^ii)n. 

(8) E E £jPj(a;iZ/t*^~^'"l'"°'''^"'modA^) 

1=0 k=0 

x\xiy^ mod N){xiy^^'^ mod N\, (g = nj-l). (51) 
One then further calculates the contribution of the operator term (51) to the 
total intensity I{m) of the multiple-quantum spectrum of the subensemble 
Sn of the density operator p(m) (49). When m ^ k'r^i, {k' — 0,1,...;/ = 
0, 1, t — 1), the contribution is given by 

t-lT-xj-l N-l 

E E ||p,-(W-''"^'"°'''"^'niod7V)P = |£, E Pj{kf. 

1=1 k=0 k=l 

where the diagonal operator term Pqq = |0)(0| is not included. When m = k'r 
the contribution is nothing to the multiple-quantum spectrum, that is, the 
integers m = k'r are the zero points of the multiple-quantum spectrum. 
When m = k'r^^ but m ^ k'r the contribution of the term (51) satisfies 

'^^iH E ||Pj-(xzy['=-«'"]'^°^'"==*modA^)|2 

1=1 k=0 

N-l 

>HE + 

A;=0 

where \pj{k)\^s.x is the maximum diagonal element of the initial density 
operator component Ps„{^)j- Now the total intensity /(m) of the multiple- 
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quantum spectrum can be calculated from the density operator (49) using 
the orthogonal relations between any pair of operator terms including the 
operator term (51) in the density operator p(m) (49). When m ^ k'vxi, the 
total intensity I{m) is given by 

ni N-l 

i{m) = i2«-i£, y: e pM- 

j=l k=l 

When m = k'r the total intensity is zero, that is, 

J(m) = if m = 0, r, 2r, 
indicating that the integers m = k'r are the zero points of the multiple- 
quantum spectrum. When m — k'r^i but m ^ k'r the total intensity I{m) 
satisfies 

7(m) > i2"-i£, E { E Pi{kf -{p + q- l)\pj{kYU.}. 

j=l k=0 

Obviously, the initial density operator component pg {0)j of the initial den- 
sity operator p(0) (47) can be efficiently transferred into the multiple-quantum 
coherences under the conditional unitary operation Ui^.s„{y, r, N)"^ with any 
integer m except the zero points m — k'r when the density operator compo- 
nent Ps„{0)j satisfies, 

» {p + q-i)\Pj{kyUs.. (52a) 

k=0 

and 

7V-1 L-1 

{ E Pjikr}-'{ E Pjik)'} - polyin) (526) 

Using the auxiliary experiment p{—m) = U/^ ^^{y, r, A^)™p(0)f//^.5^(y, r, N)"^ 

one can separate the antisymmetric part Py{m) {Py{m) = ^p{m) — |p(— m)) 
from the density operator p{m) (49). The density operator Py{m) is a pure 
multiple-quantum coherence operator. For example, the operator term (51) 
contains a component of the density operator Py[m): 

Plim) = K^ii)i(8)-(8)(^ii)i-i(8)^ix(8)(£^ii),+i(8)-0(^ii)n, 

E E £j[Pj(x,y['=-«'"l'^°'^''-'modA^) -pj(x,y[*^+«"*+'"l'^°'^''-^modA^)] 

1=0 k=0 

x{\xiy'^'^"^ mod N){xiy'' mod N\ + \xiy'' mod N){xiy'^'^^ mod N\} 
+3(^11)1 (8) - <S}{Eii)j-i (g) Ijy <S}{Eii)j+i (g) ... (8)(^ll)n, 

t~irxi-l 

E E ej[pj{xiy^''-'i"'^'^°'^''-imodN) 

1=0 k=0 

+p^.(a;,y['^+9™+™l mod N)]{i\xiy^ mod N) {xiy''+"' mod N\ 

-i\xiy''~^'^ mod N){xiy'' mod N\}, {q = n^-l). (53) 
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One sees that the zero points of the density operator Py{m) are m = k'r/2, 
that is, Py{k'r/2) = 0. When m ^ k'r,j.j2^ it is easy to calculate the contri- 
bution of the operator term Py{m) (53) to the total intensity Iy{m) of the 
multiple-quantum spectrum of the density operator Py{m), 

N-l 

mm)^TT{\pl{m)\''}^\e, E pA^?- 

fe=i 

When m — k'rxj2 but m ^k'r/2 the contribution of the operator term (53) 

N-l 

satisfies, I^{m) > 5^ p^(A;)^ — (p-|-g — l)|p^(A;)^|inax}- The total intensity 

k=l 

Iy{m) then can be calculated through the density operator Py{m). It is easy 
to prove that the intensity Iy{m) satisfies, 
Iy{m) = if m = 0, r/2, r, 3r/2, 

iV-l 

Iy{m) > |£i2-'-HE Pj{k)'-{p+q-l)\pAkfU..}iimj^0,r/2,r,3r/2,.... 

k=0 

Therefore, the pure multiple-quantum coherences of the density operator 
Py{m) are efficiently created by the conditional modular exponential opera- 
tion [^/„^5„(y, T, A^)™ with any integer m except the zero points m = k'r/2 
when the initial density operator component pg^{0)j satisfies the conditions 
(52a) and (52b). 

In order to calculate analytically the multiple-quantum spectrum of the 
density operators p(m) and Py{m) a compact and analytical derivation for 
the density operators is given below. With the help of the Fourier transform 
(12b) and the Hermitian property of the density operator the initial density 
operator p(0) (47) can be written in a symmetrical form 

m t — l^'x; — It's; — 1 

j=l 1=0 s=0 s'=0 

ni L— 1 

+Pj{s'-s,xi)\^A^i)){'^s{xi)\}+ E ejIj.<S> E PjikmC'l (54) 

j=l k=N 

where the coefficient Pj{s — s', xi) is given by 

pj{s - s' , xi) = pj{xiy''modN)exp[-i2TTk{s - s')/r^^]. (55) 

Now the time evolution of the spin ensemble is calculated with the aid of 
the Hamiltonian of Eq.(15) and the eigen-equation (4a) when acting the 
conditional unitary operation Ui^,s„{y,r, N)^ on the density operator p(0) 
ofEq.(54), 

P{m) = C/7„,5„(y,r,7V)X0)C/,+ 5jy,r,7V)- 
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= i E E E exp[-i7inim{s - s')/r^J exp[i27rm(s - s')/r^J^] 

1=0 s=0 s'=0 
ni 

x{E - s',xi){IjxCOs[27rms'/rxi\ - IjySm[27rms' /rx^\)} 

(g) + I E E E expH7rnim(s' - s)/r^J 

i=0 s=0 s'=0 

x{E ^iPj(s' - s,a;i)(Jj;cCos[27rms7'^a;J - IjySm[2TTms' /r^i])} 
X exp[i27rm(s' - s)/r^/^] (g) 

ni L— 1 

+ E^.^.x(8) E (56) 

j=l k=N 

The unitary diagonal operator exp[i27rm(s — s')/rxilz] can be expanded in 
the LOMSO subspace [29] of the subensemble 

exp[i27rm(s - s')/r^/J = ao-^o + "i-P^i + "2^2 + ••• + an.K, (57) 
where the operator Fk is the fuU symmetrical /c— body interaction basis op- 
erator of the LOMSO subspace [21, 29], 

rii fii 

Fq = E, Fi = Iz, -^2 = E '^hzhz-, F3 = E ^hzhzlmz, 

l>k=l m>l>k=l 

F = 2"^~^/'i To T 
and the coefficient ctp is generally expressed as 

rii 

Q^P^E Cpgexp[i27rm(s-s')/r^,(ni/2-g)], (58) 

9=0 

where the real coefficient Cpq can be determined using the method in Ref.[29]. 

By inserting Eqs.(55), (57), and (58) into Eq.(56) and then using the Fourier 
transform (12a) and dividing the LOMSO operator Fp into two parts: Fp = 
Flp + 2IjzF2p the density operator p{m) can be expressed in terms of the 
conventional computational basis, 

t—l m ni 

p(m) = -i X) E CpqejIjyQi{j,p,xi,k,7n) 

1=0 k=0 p,q=Oj=l 

<8)(^){+k«y^+''"' mod N) {xij/+'"''+"' mod iV| 
- mod N) mod N\} 

t—l — 1 m ni 

+1 E E ^ II CpqejIj^Qi{j,p,xi,k,m) 

1=0 k=0 p,q=Oj=l 

(g)(|){+|x/7/'=+9"^ mod N){xiy^+'i'^+'^ mod N\ 
^^^yk+qm+m j^q^j {xiy'^+i"' mod N\} 



26 



ni L—1 

+ EejI,^<S> E Pjik)\k){k\, (59) 

j=l k=N 

where the operator function Qi{j,p,xi, k,m) is defined as 
Qi{j,P,xi,k,m) = Fip[pj{xiy'' mod N) + p^(xi?/^+™modiV)] 
+F^p[pj {xiy^ mod N) - pj {xiy^~^^ mod N)] 
Using the auxihary experiment p(-m) = U^„.s„{y, N)"'p{0)Ui^.s^{y, r, N)"" 
one can further separate the antisymmetric part Py{m) from the density op- 
erator p(m), {Py{m) = \p{m) — \p{—m)), 

Py{m) = -i E E EE CpgejIjyQi{j,p,xi,k,m) 

1=0 k=0 p,q=Oj=l 

^{^){+\xiy''+^"'modN){xiy^+i"'+"' mod A^l 
-\xiy''+i"'+"' mod N) {xiy'^+'i'^ mod N\} 

t—l m ni 

-i E E ^ ^ Cpq^jIjyQ2{j,P,xi,k + qm,m) 

1=0 k=0 p,q=Oj=l 

^{^){+\xiy'' mod N) {xiy^^"^ mod N\ 
— \xiy^^'^ mod N) {xiy^ mod N\} 
t—l ^xi—l m m 
+1 E E 11 Cpq£jIjxQi(j,P,xi,k,m) 

1=0 k=0 p,q=Oj=l 

(g)(|){+|a;i|/*=+''"' mod N) {xiy''+i"'+"' mod N\ 
^^^^yk+qm+m {xiy'^+i"' mod N\} 

t—l ^xi^l m rii 

-i E E CpgSjIj^Q2{j,p,xi,k + qm,m) 

1=0 k=0 p,q=Oj=l 

mod A^) mod N\ 

+ \xiy''+'^'' mod N){xiy'' mod N\} (60) 
where the operator function Q2{j,p,Xi, k^m) is defined as 
Q2{j,p,xi,k,m) = F(p[p J {xiy'' mod N) + pj{xiy''^"^ mod N)] 
-Fip [p^ {xiy'' mod N) - pj {xiy^^^^ mod N)] 
Obviously, the density operator Py{m) is antisymmetric, that is, Py{kr/2 + 
m) = —py{kr/2 — m) {k = 0, 1, ...), and it is a pure multiple-quantum co- 
herence operator of the subensemble Sn- The antisymmetric property might 
be helpful for speeding up the searching for the zero points of the density 
operator p„(m). 



5. Searching for the period of modular exponential function 
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In previous sections it has been shown that the modular exponential op- 
eration can be performed easily on an NMR quantum computer just like 
on a classical digital computer. The classical computer outputs the value 
of the modular exponential function f{y,m,N) = y"^modN given the in- 
put integers y, m, and N, while the NMR quantum computer outputs the 
multiple-quantum spectrum of the spin ensemble which intensity does not 
reduce exponentially as the qubit number of the spin ensemble. Both the 
classical and quantum computations of the modular exponential function 
have the same computational complexity. In classical computation the val- 
ues of the modular exponential function with different inputs, i.e., the integer 
m, are generally independent on each other. Therefore it is a hard problem 
to find the period r of the modular exponential function on a classical com- 
puter. However, the essential difference for the quantum computer from the 
classical one is that the computational process on the quantum computer 
obeys the unitary dynamics of quantum mechanics. Then in the factoring 
problem the quantum computational process and output (through the den- 
sity operator) are governed by the Liouville-von Neumann equation or the 
Schrodinger equation where the integer m acts as the discrete time variable, 
as can be seen below, and therefore the output results at different times (m) 
really correlate to each other. This essential point could form the base to 
solve efficiently the factoring problem and play a key important role for the 
quantum computer outperforming the classical one in solving the factoring 
problem. 

In general, the amplitudes and phases of multiple-quantum coherences of 
p(m) and Py{m) with different quantum orders are dependent on the inte- 
ger m in the spin ensemble, which are described by the Liouville-von Neu- 
mann equation. In particular, the amplitude and phase for the long-range- 
interaction and higher-order multiple-quantum coherences could be helpful 
for efficiently searching for the period r. The Liouville-von Neumann equa- 
tion with Hamiltonian H{y, r, N) that governs the unitary dynamical process 
during the modular exponential operation (1) in a spin ensemble can be writ- 
ten as 

dp{t)/dt = -t[H{y,r,N),p{t)], (/i= 1). (61) 
The solution to the Liouville equation then is given formally by 

p{t) = U{t)p{0)U{t)+ (62) 
where the propagator is written as 

U{t) = exp[-itH{y, r, N)]. (63) 
By comparing the unitary operator U{t) with the modular exponential op- 
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erator U{y,r, N)"^ of Eq.(7) one sees that the integer m in the modular 
exponential operator is really equivalent to the time variable t and their dif- 
ference is merely that the time variable t is continuous but the integer m 
discrete. Therefore, the propagator is time periodic: U{t) = U{t + r), where 
the period r needs to be determined in the factoring problem. Below it is 
assumed that the Hamiltonian H{y,r,N) consists of a dominating and a 
relative small operator components, 

H{y,r,N) = Ho + Hi. (64) 
The dominating term Hq is a specific order quantum operator, for exam- 
ple, a zero-quantum coherence operator, while the small operator term Hi 
is usually a multiple-quantum operator. By making the coordinate frame 
transformation: p^(t) = exp{iHQt)p(t) exp^—iHot), here the frame is called 
the interaction frame defined by the Hamiltonian Hq, the Liouville equation 
(61) is rewritten as 

dp,{t)/dt^-i[Hi{t),p,{t)], (65) 
with the time-dependent Hamiltonian in the interaction frame: 

Hi{t) = exp{iHot)Hi exp{-iHot). (66) 
The solution to the Liouville equation (66) is given in form 

p,{t) = Ui{t,to)p{to)Ui{t,toy (67) 

with the propagator in the interaction frame: 

t 

Ui{t,to) =Texp{-i J Hi{t')dt'), (68) 

to 

where the operator T is Dyson time-ordering operator. To see more clearly 
the time evolution process in the interaction frame the solution of Eq.(67) is 
expanded, 

p{t)^exp{-iHot){p{0)-it[Hi{t),p{0)] 

-l^lHiit), [Hi{t),Pm] + -ycxpizHot). (69) 
Suppose that the initial density operator p(0) is a LOMSO operator, for 
example, the initial density operator component Ps„{^) of Eq.(46) that is a 
LOMSO operator of the subensemble Sn- It can be seen clearly from the 
expansion (69) how the initial density operator is converted into multiple- 
quantum coherences as time development. Since the zero-quantum Hamilto- 
nian Hq is dominating the initial LOMSO density operator is converted effi- 
ciently into the zero-quantum coherence at a short time, e.g., t = 1 (m = 1). 
In general, the p— order quantum peak is strongest at a short time (a small 
m) if the Hamiltonian H{y,r,N) contains a dominating p— order quantum 
coherence operator. This ensures that the NMR multiple-quantum signals 



29 



at the time points t — and t — 1 can be precisely distinguished experimen- 
tally without an exponential resource by the multiple-quantum spectroscopic 
method. The measurement precision is important on a quantum computer, 
while it is not any problem in a classical computer. Since the multiple- 
quantum Hamiltonian Hi{t) is small the nonzero-order multiple-quantum 
coherences will grow slowly and monotonously in a long time interval. At 
the same time the zero-quantum coherence first increases quickly and reaches 
its maximum and then decreases gradually as time development because part 
of the initial density operator is converted into the multiple-quantum coher- 
ences. The same time evolution behavior of the density operator p{t) near 
the zero point p(0) also occurs at other zero points p{kr) {k — 0, 1, ...) due to 
the period of the density operator, p{kr) = p(0). If such time development 
behavior for the zero-quantum and nonzero-order quantum coherences con- 
tinues in a time interval AT satisfying r/AT ~ poly{n) then the searching 
for the zero points of the multiple-quantum spectra will be polynomial-time 
on the NMR quantum computer. The searching efficiency is proportional to 
the time interval AT, that is. the longer the time interval AT the higher the 
efficiency. The searching efficiency will decrease if the density operator p{t) 
arrives at its steady state at a shorter time. Here the steady state implies that 
the intensity /^(p, t) for any p— order quantum coherence does not change as 
the time. The situation may occur when the multiple-quantum Hamiltonian 
Hi is not small, however, even in this case the searching is still efficient if the 
initial integer m is sufficiently near the zero points. Obviously, the searching 
for the zero points is locally efficient in a small region near the zero points. 
The steady-state problem is harmful for the present factoring algorithm to 
find efficiently the zero points in a spin ensemble and needs to be overcome. 
It is closely related to the Hamiltonian H{ij, r, N) and the distribution of 
different order quantum transitions in a spin ensemble. It is possible to over- 
come the steady-state problem by manipulating the Hamiltonian H[y,r,N) 
and choosing the proper initial density operator in the factoring sequence. 

The distribution of different order quantum transitions for a spin ensem- 
ble with n non-equivalent spins-1/2 has been found [19]. The number of 

the zero-quantum transitions is Zq — |{( ) — 2"} and for p— order quan- 

2n 

tum transitions Z„ — ( ), p = 1, 2, n. In the case of the large n 

^ n — p 

and relative small p ^ the p— order quantum transition number can be 
approximated by Stirling formulae, Zp — 4'*(7rn)~^/^ exp(— p^/n). This indi- 
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cates that the population distribution of different order quantum transitions 
is extremely nonuniform in the spin ensemble. In general, the lower order 
quantum transitions such as zero-, single-, and double-quantum transitions 
are much more than those higher-order quantum transitions in a spin ensem- 
ble. The spectral intensity of the p— order quantum transition is generally 
proportional to the probability of the p— order quantum transition in the 
distribution, although this is not absolute. Therefore, it is better to choose 
lower order quantum transition spectral peaks such as zero-, single- or double- 
quantum transitions to help the searching for the period r in the factoring 
algorithm. 

It has been shown in previous sections that the total conversion efficiency 
Iy{m)/IiQ of the multiple-quantum coherences is almost independent of the 
qubit number n when m is not a zero point and both the conditions (52a) 
and (52b) are met. Because there are only (2n -|- 1) spectral peaks in the 
multiple-quantum spectrum, on average, each peak intensity is approximately 
inversely proportional to the qubit number n even when the integer m — 1. 
Then there are at least some peaks among the 2n -f- 1 peaks, for example, the 
p— order multiple-quantum peak, which intensity Iy{p,m) can be detected 
precisely without an exponential resource. Consequently, with the factoring 
sequence in previous sections one can fix experimentally the zero points mo = 
kr/2 {k — 0, 1, ...,) from a small neighbor region (mo ± 1 at least) of the 
zero points without an exponential resource, and thus the factoring sequence 
is locally efficient in a small neighbor region of the zero points at least. 
However, it is not clear whether the factoring sequence is yet efficient or 
not when the searching for the zero points starts at those time points m 
far from the zero points. The time development behavior of the long-range- 
interaction and higher-order multiple-quantum spectral intensities may play 
an important role to find efficiently the zero points when the searching for 
the zero points starts at those points far from the zero points. The time 
development behavior is dependent on the Hamiltonian H{y, r, N) of the 
modular exponential operation. The numerical simulation using the density 
operator Pyijn) of Eq.(27) shows that the zero-quantum peak Iy{p = 0, m = 
1) is stronger than any other multiple-quantum peaks when the integer y = 
2,4, and is much smaller than N. The simulation also shows that some 
nonzero-order quantum peaks are also quite strong even for a large integer 
N and the smallest integer y — 2, although the zero-quantum peak is still 
strongest. This implies that the Hamiltonian H(y, r, A^) with a small integer 
y — 2,4, ... and a large integer N still have a quite large multiple-quantum 
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coherence component Hi in addition to the strongest zero-quantum coherence 
operator Hq, so that the strong zero-quantum peak falls off rapidly as the 
integer m and the density operator Pyim) approaches rapidly to its the steady 
state. However, the searching for the zero points using the zero-quantum 
peak or other multiple-quantum peaks is still locally efficient, that is, if the 
initial integer m is sufficiently near the zero points then the zero points can 
be found efficiently even for a large number with the help of the time 
development behavior of the zero-quantum peak or other multiple-quantum 
peaks. The searching for the zero points based on the time development 
behavior of the p— order quantum peak may become really inefficient when 
the searching starts from those points m far from the zero points. For an 
integer y ^ 2,4,..., or for a large integer y the density operator Py{m) of 
Eq.(25) approaches quickly to its steady state as the integer m. One of the 
reasons for it could be that the initial density operator p(0) of Eq.(21) is very 
special in the factoring sequence, that is, Ps„(0) is the unity operator. It can 
be known from Eq.(56) that if the initial density operator P5^(0) is a LOMSO 
operator then the multiple-quantum coherences created by the conditional 
modular exponential operation t, N)"^ belong to each isolated subset 

S{xi) X S{xi) and the maximum number (~ i(r^ — r)) of multiple-quantum 
transitions in the subset ^'(a;;) x S{xi) usually is much less than the maximum 
number |(4" — 2") of multiple-quantum transitions of the subensemble Sn 
with n spins- 1/2. In particular, the maximum number of multiple-quantum 
transitions induced by the conditional modular exponential operation on the 
initial density operator p(0) of Eq.(21) is not more than {N — 1) (|2" < 
N < 2") for any given integer m. This can be seen from the density operator 
Py{m) of Eq.(27). This number is greatly less than the maximum number 
1(4"' — 2"). This could be one of the reasons why the density operator Pyijn) 
(27) enters into its steady state rapidly as the integer m. In the future it 
will be studied in detailed how the initial density operator P5„(0) is chosen 
properly to overcome the steady state problem. 

A possible scheme to overcome the steady-state problem is described 
below from the point of view of manipulating the Hamiltonian H[y, r, N) 
of the modular exponential operation. First one finds a unitary operator 
G{y. r, N) so that the transformed Hamiltonian H{y, r, N) has a dominating 
zero-quantum coherence component: 

H{y, r, TV) = G{y, r, N)+H{y, r, N)G{y, r, N). 
Then this new Hamiltonian is acted on the initial density operator by replac- 
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ing the original Hamiltonian H(y, r, N) of the modular exponential operation, 
and the generated multiple-quantum spectra could be able to be used to effi- 
ciently find the period r. The unitary operator G{y,r,N) always exists, but 
it is a challenge how to find the exact unitary operation G{y, r, N) that can 
be implemented in polynomial time. The Hamiltonian H(y, r, N) always can 
be diagonalized unitarily. Assume that there is a unitary operator V{y,r, N) 
to diagonalize the unitary operator U{y, r, N): 

V{y, r, N)+U{y, r, N)V{y, r, N) = A{y, r, N). (70) 
The unitary operator V{y,r,N) can be constructed from the Fourier trans- 
forms (12a) and (12b), but it contains the period r in an explicit form and this 
makes it difficult to construct its explicit quantum circuit. But the unitary 
operator V{y,r,N) could be built up approximately. The Fourier transform 
(12a) over the period r^i may be replaced with the following approximated 
Fourier transform over the whole range of the integer N [2, 3, 4, 6, 7]: 

N-l 

|*,(xi))«^ E exp{i27rsk/N)\xiy''modN). (71) 

k=0 

If the period r divides the integer then the Fourier transform (71) is exact. 
But the period r usually does not divide the integer A^. Therefore, the uni- 
tary operator V{y,N) built up with the Fourier transform (71) diagonalizes 
approximately the unitary operator U {y, r, N), that is, V{y, N) « V{y, r, N). 
The unitary operator V{y, N) does not explicitly depend on the period r and 
has a polynomial quantum circuit since the Fourier transform (71) can be 
constructed efficiently [2, 3, 4, 6, 7]. By using the unitary operator V{y, N) to 
diagonalize approximately the Hamiltonian H[y,r,N) the diagonal unitary 
operator is obtained: V{y, N)^U {y"^, r, N)V{y, N). Then one chooses further 
a proper unitary operator W{y,N), which Hamiltonian has a dominating 
zero-quantum coherence component and a relative small multiple-quantum 
component, to construct the desired unitary operator: 

U{y, r, AT)- = W{y, N^Viy, NyU{y"\ r, N)V{y, N)W{y, N). (72) 
Therefore, the unitary operator G{y,r,N) that converts the Hamiltonian 
H{y, r, A^) into the desired Hamiltonian H{y, r. A") which has a dominating 
zero-quantum coherence component may be approximated by the unitary 
operator G{y, N) which does not explicitly depend on the period r and is 
given by 

G{y, r, N) « G{y, N) = V{y, N)W{y, N). ^ (73) 
Obviously, the unitary operator U{y,r,N) satisfies U{y,r, N)^ — E when 
m — kr {k — 0, 1, 2, ...,), indicating that the unitary operator tJ{y, r, N) has 
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all the periods of the original unitary operator U{y.r,N). Generally the 
unitary operator W{y, N) is chosen suitably so that the unitary operator 
U {y, r, A^) has not any other periods except the own periods of the unitary 
operator U{y,r,N). Now the new unitary operation (72) places the original 
unitary operation U{y,r,N) in the factoring sequence in previous sections. 
Then in the searching for the zero points the initial density operator p(0) 
of the spin ensemble (In^Sn) is converted efficiently into the zero-quantum 
coherence quickly at a short time, e.g., m = 1, under the conditional unitary 
operation UI„^s„{y,'^^, N)™- since the Hamiltonian H{y,r,N) has a dominat- 
ing zero-quantum coherence component. Therefore, the dominating zero- 
quantum coherence component of the Hamiltonian governs the time evolu- 
tion behavior of the spin ensemble at a short time and is responsible for 
precisely distinguishing the zero points from other time points in the time 
region near the zero points p{kr) without an exponential resource, while the 
relative small multiple-quantum coherence component of the Hamiltonian 
will be responsible for the efficient searching for the zero points starting from 
those time points far from the zero points. 

According to the factoring sequence in previous sections one may find 
a zero point m = r' by searching for the zero points of the density opera- 
tor Py{m), but the zero point r' could not be the minimum period r of the 
modular exponential function. Suppose the period r is an even integer as 
before. Obviously, the ratio r'/r can only take a half integer (2k + l)/2, 
an even integer 2k, or an odd integer [2k + 1), k = 0, 1,2, .... If the ratio 
r'/r is a half integer, i.e., r'/r = {2k + l)/2 then the modular exponential 
function f{y,r',N) = y^' modN — y^'^^'^+^^Z^ mod = y^^'^modN. One can 
use directly the function /(y, r', A^) to determine the non-trivial factor of the 
integer A^ if the function f{y, r', N) ^ —1. The non-trivial factor takes either 
gcd{f{y,r/2,N) - 1,N) or gcd{f{y,r/2,N) + 1, A^) [1, 2, 3]. Therefore, one 
needs merely to find a half- integer zero point r' satisfying r'/r = {2k + l)/2 to 
factor the integer A'^. Suppose that the factoring sequence finds a zero point 
r'. One uses the zero point r' to calculate the function f{y, r', A^) which will 
take about 0{{log2Nf) steps [1, 2]. If f{y,r',N) = 1 then r'/r = 2k or 
2k + 1, otherwise r'/r = {2k + l)/2. For the case r'/r = 2k ot 2k + 1 one 
further calculates the function f{y,r'/2,N) by using the integer r'/2, that 
is, f{y,r'/2,N) = i/*"'/^ mod A". If now f{y,r'/2,N) 7^ ±1 one can use the 
function f{y,r'/2,N) to find further a non-trivial factor of the integer A^; 
otherwise f{y,r'/2,N) = 1 and r'/2 — kr. For the case r'/2 = kr one calcu- 
lates the function f{y,r'/4,N) again. If f{y,r'/4,N) ^ ±1 one will obtain 
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a correct function f{y,r' /4, N) to factor the integer N, otherwise calculate 
further f{y,r'/8,N). Therefore, hj p = 0(n) steps at most to calculate 
the function /(y, r'/2'^, A^), A; = 0, 1, 2, — 1 one can finally find a correct 
function f{y,r' /2p~^, N) to factor the integer N. If the period r is not an 
even integer or f{y, r', N) — y^ mod A?" = —1, meaning that one can not find 
a non-trivial factor of N by the function f{y, r', N), one needs to choose an- 
other integer y coprimc to the integer [1, 2, 3] and then run the factoring 
sequence above to find a zero point r' so as to obtain the correct f{y, r', N). 



6. Discussion 

In this paper a quantum factoring sequence based on the unitary dynam- 
ics of quantum mechanics has been proposed to solve the prime factorization 
problem on a spin ensemble without any quantum entanglement. It uses 
the NMR multiple-quantum measurement techniques to output its quan- 
tum computational results. The NMR quantum computer can perform the 
modular exponential operation just like a classical digital computer, but its 
quantum computational output is the inphasc multiple-quantum spectrum 
of the spin ensemble which may reduce merely in a polynomial form as the 
qubit number of the spin ensemble. The computational complexity of the 
modular exponential operation is the same on both the quantum computer 
and the classical one. Quantum entanglement is not involved in the present 
ensemble quantum computation of prime factorization because there is not 
any quantum entanglement in the spin ensemble used to perform the prime 
factorization. The time evolution process of the modular exponential opera- 
tion on the quantum computer obeys the unitary dynamics of quantum me- 
chanics and hence the computational output is governed by the Liouville-von 
Neumann equation of quantum dynamics. This essential difference between 
the quantum computer and the classical one could be the key point for the 
quantum computation outperforming the classical one in the prime factor- 
ization on a spin ensemble without any quantum entanglement. It has been 
shown that the prime factorization based on the unitary dynamics of quan- 
tum mechanics on a spin ensemble is locally efficient at least. Therefore, the 
quantum entanglement could not be a unique resource to achieve speedup 
of quantum computation in the prime factorization on a spin ensemble and 
quantum dynamics could play an important role for the origin of power of 
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quantum computation. The steady-state problem is a harmful problem. It 
hampers the present factoring sequence to find efficiently the period of the 
modular exponential function. It is worth studying in detailed in the future 
how the steady-state problem is dependent on the initial density operator of 
the factoring sequence and how the steady-state problem may be overcome 
by manipulating the Hamiltonian of the modular exponential operation. 

There are a number of works [30, 31] to describe how to construct effi- 
ciently the quantum circuit of the modular exponential unitary transforma- 
tion U {y, r, N) in a quantum system with qubit number much more than 
(1 -|- [log2 N]), where a large number of extra auxiliary qubits are used. The 
construction of the unitary operator U {y, r, N) may be easier in a quan- 
tum system with a larger Hilbert space, that is, with a larger number of 
qubits. However, in practice it is still a challenge to construct efficiently the 
quantum circuit of the unitary operator U {y, r, N) in a spin ensemble with 
(1 -I- [log2 A?"]) qubits at least. The implementation for the Shor's factoring 
algorithm on a quantum system need consume a number of qubits, but if the 
modular exponential operation could be implemented efficiently on a spin 
ensemble with qubits as low as (1 + [log2 A^]) then this would simplify greatly 
the implementation of the prime factorization. 

Multiple-quantum coherences are generally measured indirectly through 
the detection of single quantum coherence in NMR spectroscopy. The mea- 
surement is more time-consuming than the direct detection of single quantum 
coherence. However, the importance is that the measurement time for each 
running of the factoring sequence based on the NMR multiple-quantum spec- 
troscopic method is almost independent of qubit number of a spin ensemble. 
If digital resolution to record experimentally NMR multiple-quantum signal 
needs to keep constant then the consuming time is approximately linearly de- 
pendent on the qubit number since the spectral width to cover over all 2n + 1 
multiple-quantum spectral peaks is about 2nu!s, approximately proportional 
to the qubit number n. Therefore, the measurement is not a severe computa- 
tional complexity problem in the factoring sequence. An improved method to 
overcome the time-consuming problem of multiple-quantum coherence indi- 
rect measurement might be using one-dimensional multiple-quantum filtering 
experiments in NMR spectroscopy [19]. The one- dimensional experiments 
should employ gradient magnetic field [32] instead of the phase cycling to 
select the multiple-quantum coherence with desired quantum order before 
detection and then convert it into single quantum coherence to be detected 
directly. In the one-dimensional experiments the antisymmetric property of 
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the density operator Py{m) might be useful for the speedup of the searching 
for the zero points of the density operator. However, there are some prob- 
lems to be solved for the one-dimensional experiments to be used in the NMR 
quantum computation, for example, how to convert efficiently the desired or- 
der quantum coherences into inphase single quantum coherence which can 
be detect efficiently. 

Relaxation or decoherence effect in a spin ensemble is usually harmful 
for any ensemble quantum computation, but it might be harmless for the 
searching for the zero points in the factoring sequence on a spin ensemble. 
Since the NMR multiple-quantum coherences usually decay in an exponential 
form and irreversibly as the time development in a spin ensemble then the 
time development behavior of the multiple- quantum spectral peaks such as 
the zero-quantum peak used to search for the zero points may become more 
distinct in the region near to the zero points so that the searching might 
become more efficient. However, the decoherence effect may destroy the 
efficient detection for multiple-quantum coherences and especially for those 
higher-order quantum coherences due to the fact that a higher-order quantum 
coherence usually has a shorter relaxation time and its NMR signal usually 
decays much more rapidly than those lower order quantum coherences in a 
spin ensemble. Therefore, relaxation effect is a compromised effect on the 
present prime factorization on an NMR quantum computer. 
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